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SUPERNOVA 1987A: A TEMPLATE TO LINK SUPERNOVAE TO THEIR REMNANTS 
S. Orlando^ M. MicelE’\ M.L. Pumo^’^’^, F. Bocchino^ 

ABSTRACT 

The emission of supernova remnants reflects the properties of both the progenitor supernovae and 
the surrounding environment. The complex morphology of the remnants, however, hampers the 
disentanglement of the two contributions. Here we aim at identifying the imprint of SN 1987A on 
the X-ray emission of its remnant and at constraining the structure of the environment surrounding 
the supernova. We performed high-resolution hydrodynamic simulations describing SN 1987A soon 
after the core-collapse and the following three-dimensional expansion of its remnant between days 
1 and 15000 after the supernova. We demonstrated that the physical model reproducing the main 
observables of SN 1987A during the first 250 days of evolution reproduces also the X-ray emission of the 
subsequent expanding remnant, thus bridging the gap between super novae and supernova remnants. 
By comparing model results with observations, we constrained the explosion energy in the range 
1.2 — 1.4 X 10^^ erg and the envelope mass in the range 15 — 17M0. We found that the shape of X-ray 
lightcurves and spectra at early epochs (< 15 years) reflects the structure of outer ejecta: our model 
reproduces the observations if the outermost ejecta have a post-explosion radial profile of density 
approximated by a power law with index a = — 8. At later epochs, the shapes of X-ray lightcurves 
and spectra reflect the density structure of the nebula around SN 1987A. This enabled us to ascertain 
the origin of the multi-thermal X-ray emission, to disentangle the imprint of the supernova on the 
remnant emission from the effects of the remnant interaction with the environment, and to constrain 
the pre-supernova structure of the nebula. 

Subject headings: hydrodynamics — instabilities — shock waves — ISM: supernova remnants — X- 
rays: ISM — supernovae: individual (SN 1987A) 


1. INTRODUCTION 

Supernova remnants (SNRs), the leftovers of super¬ 
nova (SN) explosions, are diffuse extended sources with 
a quite complex morphology. General consensus is that 
such morphology reflects, on one hand, the physical and 
chemical properties of the progenitor SN and, on the 
other hand, the early interaction of the SN blast wave 
with the inhomogeneous circumstellar medium (CSM). 
Thus investigating the intimate link that exists between 
the morphological properties of a SNR and the complex 
phases in the SN explosion may help: 1) to trace back 
the structure and chemical composition of SN ejecta, and 
the dynamics and energetics of the SN explosion; 2) to 
probe the structure and geometry of the CSM, thereby 
mapping the final stages of the stars evolution. How¬ 
ever, the very different time and space scales of SNe 
and SNRs m ake difficult to study in detail their con¬ 
nection fe.g. [Badenes et al.l l2QQ8l: lEllinger et al.l I2Q12I: 
lYamaguchi et al.ll2Ql4 iPatnaude et al.ll2Q15[)~ 

Because of its youth and proximity, SN 1987A in the 
Large Magellanic Cloud offers a unique opportunity to 
bridge the gap between the progenitor SN and its rem¬ 
nant. SN 1987A was an hydrogen-rich core-collapse 
(CC) SN tha t was observed in outburst on Eebru- 
ary 23, 1987 (|West et al.l Il987l) . It occurred approx- 
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imately 51.4 kpc from Earth (|Panagial Il999[) . About 
80 days after the explosion, it reached a peak appar- 
ent visual magnitude of ^ 3 (naked eye visible; e.g. 
ICatchpole et al.l 119881: iHamuv et al.l Il988[) . Its evolu¬ 
tion has been accurately monitored in different wave¬ 
length b ands since the out burst, from infrared (e.g. with 
Spitzer^ iDwek et al.l |2Q1Q11. to optical (e . g. with Hub¬ 
ble S pace Telescove, iLawrence et al.|[2QQQl: iLarsson et al. 
2011 ). to X-ray bands (e.g. with Rosat ARdheil etTT 
2006. Chandra, [Helder et M] 120131 and XMM-Newton, 
Haberl et al.ll2006l) . This has provided a wealth of high- 
quality data with unprecedented completeness, making 
of SN 1987A an ideal template to study the SN-SNR con¬ 
nection. 

The study, however, is complicated by the interac¬ 
tion of the blast wave with the surrounding inhomo¬ 
geneous medium. Optical images soon after the out¬ 
burst re vealed an enigma tic triple-ring nebula around 
the SN (|Crotts et al.lll98^ . The nebula consists mainly 
of a dense central equatorial ring and two outer rings 
displaced by about 0.4 pc above and below the cen¬ 
tral ring and lying in planes almost parallel to the 
equatorial one. It has been suggested that the nebula 
might be the result either of the merging of two massive 
stars occurred some 20,000 ye ars before the explosion 


(|Morris & Podsiac 

lowskil 12007 

) or of mass loss from a 

fast-rotating star ( 

Chita et al.l 

200l. 


The interaction of the SN with the nebula is best ob¬ 
served in the radio and X-ray bands. It started about 
3 years after the explosion when both radio and X-ray 
fluxes began to incre ase with time (|Hasinger et al .111 9961: 
iGaensler et al.l Il997[ ) . This was interpreted as due to 
outer ejecta lighting up the dense wind emitted during a 
previous red supergiant (RSG) phase and subsequently 











































2 


S. Orlando et al. 


swept-up by the fast wind during a phase of blue super¬ 
giant (BSG) (|Chevalier fc Dwarkadaslll995[) . After about 
16 years, suddenly the soft X-ray light curve has steep¬ 
ened still fu rther, contrary to the hard X-ray and radio 
lightcurves (jPark et al.ll2QQ5[ ). This was interpreted as 
due to the blast wave swee ping up the dense central equa¬ 
torial ring (jMcCravll2QQ7[ ) . Since then the soft X-ray flux 
continued to grow, indicating that currently the shock is 
still traveling through t he dense part of the equatorial 
ring (|Helder et alJl2Ql^ . 

X-ray observations more than others encode informa¬ 
tion about the physical properties of both the nebula and 
the stellar ejecta. Decoding these observations, there¬ 
fore, might open the possibility to reconstruct the neb¬ 
ula and ejecta structures soon after the SN explosion. 
This requires accurate and detailed numerical models 
able to follow the ejecta evolution from the SN explo¬ 
sion to the SNR phase and to reproduce altogheter the 
emission properties of both the progenitor SN and its 
remnant. 

Current models, however, were us ually aimed at de¬ 
scrib i ng either the SN e volution (e.g. iPumo fc Zampieril 
l2Qlll: lHandv et aL[l2Q14 ) or the expansion o f the remnant 
(e.g. IZhekov et al1l2QQ9 : iDewev et al.ll2Ql3 iPotter et al.l 
12014 ). The former models describe the early SN evolu¬ 
tion and ignore its subsequent interaction with the neb¬ 
ula; the latter assume an initial parametrized ejecta pro¬ 
file that is not proven to reproduce the main observables 
of the progenitor SN, leaving out an accurate description 
of the ejecta distribution of mass and energy soon after 
the SN explosion. These limits make difficult to disen¬ 
tangle the effects of the initial conditions (i.e. the SN 
event) from those of the boundary conditions (i.e. the 
interaction with the environment). A first attempt to 
connect some properties of CC-SN (e.g. the composition 
and a parameterization of the circumstellar environment) 
to some observa ble bulk properties of S NR has been dis¬ 
cussed recently (|Patnande et al.ll2Q15f ). even though by 
adopting a ID approach. 

The structure of the pre-SN nebula of SN 1987A is 
inherently three-dimensional (3D), as is visible in the 
images from the H ubble Space Telescope (HST) (e.g. 
iLarsson et al.lDOlU ) . A proper description of the interac¬ 
tion of the ejecta with the nebula therefore requires 3D 
calculations. To da te only a 3D model has been proposed 
(|Potter et al.ll2Ql4) to explore the origin of the asymmet¬ 
ric radio morphology observed in SN 1987A. However this 
model focusses on the analysis of the radio emission and 
its initial condition is a parametrized function describ¬ 
ing the ejecta distribution about two years after the ex¬ 
plosion. As a consequence their model cannot describe 
the evolution of the SN and their adopted distribution of 
ejecta is not proven to be consistent with the observables 
of the SN. 

Here we present a hydrodynamic model describing the 
evolution of SN 1987A from the breakout of the shock 
wave at the stellar surface until the expansion of its 
remnant through the nebula surrounding the SN. We 
ran several high-resolution simulations to reproduce al¬ 
together the main observables of the SN (i.e. bolomet- 
ric lightcurve, evolution of line velocities, and continuum 
temperature at the photosphere) and the properties of 
X-ray emission of the remnant (lightcurves, spectra, and 
morphology). Our simulations cover the first 40 years 


of evolution to make predictions on the future evolution 
of the remnant in view of the spatially-resolved high- 
resolution spectroscopy capability of the for thcoming X- 
ray observatory Athena (|Nandra et al.l[2QT^ . The paper 
is organized as follows. In Section [2] we describe the hy¬ 
drodynamic model and the numerical setup, in Section [3] 
we describe the results and, finally, we draw our conclu¬ 
sions in Section m 

2. PROBLEM DESCRIPTION AND NUMERICAL SETUP 

We first modeled the SN by performing one¬ 
dimensional (ID) simulations of the “early” post¬ 
explosion evolution of a CC-SN during the first 250 days 
(see Section [2T]) . Then the output of these simulations 
was mapped into 3D and provided the initial condition 
for the ejecta structure of a 3D model describing the fur¬ 
ther evolution of the SN and the subsequent development 
of the SNR interacting with the CSM (see Section [2^ . 

2.1. Modeling the post-explosion evolution of the 
supernova 

The post-explosion evolution of the SN was mod¬ 
eled by using a ID Lagrangian code, specifically tai¬ 
lored to simulate the main observables in CC-SNe 
(namely the bolometric lightcurve and the time evo- 
lution of the pho t osphe ric velocity and temperature; 
iPumo fc Zampier IMnL an d wid e ly used to model 
observed events dSpiro et al.l l2Ql4 iTakats et al.l l2Ql4 
iDairOra et al.l 12014) i ncluding the SN 1987A-hke ones 
(|Pastorello et al.ll2Q12[) . 

The code solves the equations of relativistic radia¬ 
tion hydrodynamics in spherical symmetry for a self- 
gravitating matter fluid which interacts with radiation. 
Its distinctive features are: 1) a fully general relativistic 
treatment; 2) an accurate treatment of radiative trans¬ 
fer at all regimes (from the one in which the ejecta are 
optically thick up to when they are optically thin); 3) 
the coupling of the radiation moment equations with 
the equations of relativistic hydrodynamics during all 
the post-explosive phases, adopting a fully implicit La¬ 
grangian finite difference scheme; 4) the computation of 
heating effects due to decays of radioactive isotopes syn¬ 
thesized in the SN explosion; and 5) the computation of 
the gravitational effects of the central compact object on 
the evolution of the ejecta. 

Thanks to these characteristics, the code is able to 
compute the evolution of the stellar ejecta and emitted 
luminosity from the breakout of the shock wave at the 
stellar surface up to the so-called nebular stage (i.e. when 
the envelope has recombined and the energy budget is 
dominated by the radioactive decays of the heavy ele¬ 
ments synthesized in the explosion, see also Section [3T|) . 
At the same time, the code is able to accurately deter¬ 
mine the fallback of material on the central object and, as 
a consequence, the amount of ^^Ni present in the ejected 
envelope at late times. 

The simulations start from the breakout of the shock 
wave at the stellar surface, using initial conditions 
that mimic the physical properties of the ejected ma- 
terial after shock passa ge following the core-collapse 
(jPumo fc Zampie^l2Qll[ ). As a consequence, the mod¬ 
els depend on some basic parameters that characterize 
the main radiative, thermal, and dynamical properties of 
such material. These parameters are: the progenitor ra- 
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dius Rq, the total ejecta energy L’sn, the envelope mass 
at shock breakout Menv, and the total amount of ^^Ni 
initially present in the ejected envelope Mni- In our sim¬ 
ulations Menv coincides with the total mass surrounding 
the central compact object at the start of simulations. 
So, during the post-explosive evolution, Menv is equal to 
the ejected mass Mej plus the mass fallen back to the 
central compact object. In our models of SN 1987A, the 
latter mass is of the order of a few hundredths of a solar 
mass and, as a consequence, Menv essentially corresponds 
to Mej. 

We performed several simulations of the SN evolution, 
exploring the parameter space defined by Menv and L^sn • 
The exploration has been restricted to the ranges of val- 
nes commonly discussed in the literature for SN 1987A 
(lArnettI 11987c IWooslev 19881: lUtrobin fc Chngail I2QQ5I : 
iPnmo fc ZamDierill2Qlll: iHandv et al.ll2Q14l) : we consid- 
ered models with Menv ranging between 8 and 19 Mq 
( including the bulk of ejecta and the high-velocity shell), 
and Esn ranging between 1 and 2 x 10 ^^ erg. Also, in 
our model, the outermost high-veloc ity shell (4000 < v < 
20000 km s“^ ; iLawrence et al.ll 2000 f) can be described by 
a power law with ind ex a (|Shigevama fc Nomotol [19901: 
iBlinnikov et al.l[20001 ). Here we explored models with a 
ranging between —8 and —9, a ccording to the values sug¬ 
gested in the literature (e.g. IWooslevlll988[) . We fixed 
all the other parameters of the model, namely the initial 
radius Rq = 3 x 10^^ cm (w hich is a reliable value for 
the progenitor of SN 19 87A; lArnettI [19961: lYonnd [2004 
iPnmo fc ZamDierill201l[) and the initial amount of ^^Ni 
Matz = O.O 7 M 0 (corresponding to the amo unt of Ni syn¬ 
thesi zed during the explosion of SN 1987A; I Arnett et al.l 

[HU). 

Note that our ID SN model cannot simulate possible 
asymmetries developed during the explosion as suggested 
by a significant body of observational and simulation evi¬ 
dence. It goes without saying that the most desirable way 
to model a CC-SN would include a multi-dimensional 
hydrodynamical calculation, following the core-collapse, 
the bounce at nuclear densities, the leakage of neutri¬ 
nos from the proto-neutron star, the neutrino heating, 
and the delayed explosion with the shock propagating 
through the envelope to be ejected. However, current 
multi-dimensional models can follow the evolution of the 
exploding star only for a very short time^ (of the order 
of few minutes), preventing any possibility to use their 
output for post-explosion calculations that follow in de¬ 
tail the evolution of the CC-SN ejecta. Conversely, our 
ID model allows us to follow the post-explosion evolu¬ 
tion of the SN for days, to describe in detail the fallback 
by using a full-relativistic approach and, therefore, to 
estimate accurately the total mass ejected in the SN ex¬ 
plosion. Also our lagrangian code is not flux-limited (as 
the multi-dimensional codes are), allowing us to simulate 
accurately the evolution of the ejected material during 
both the initial phase of the shock breakout and the fol¬ 
lowing nebular phase in which the ejected SN envelope 
has recombined. 

2 . 2 . Modeling the supernova remnant evolution 


^ However note that a new approach to follow the post-explosion 
evolution of the SN i n multi-dimensions fo r a longer time has been 
proposed recently bv IJoggerst et aP (|2009l b 


Optical images clearly show that the structure of the 
nebula as well as the m orphology and evolution of the 
ejecta are in herently 3D (jKjaer et al.lDOlOl: iLarsson et al.l 
l2QllL [201^ . Therefore, once the early phase of the SN 
explosion has been modeled in ID (see Section 12 . 1 [) . the 
output of the SN simulation was mapped into 3D. Then 
the subsequent evolution and transition from the SN 
phase to the SNR phase were modeled by numerically 
solving the time-dependent fluid equations of mass, mo¬ 
mentum, and energy conservation; the equations in a 3D 
Cartesian coordinate system (x, z) take into account 
the radiative losses from an optically thin plasma: 


dt 


V • pu = 0 , 


dpu 

dt 


+ V • puu + VP = 0 


dpE 

dt 


+ V • {pE -h P)u = -n^nuEiT) , 


( 1 ) 


where E = e-h |up/2 is the total gas energy (internal 
energy, e, and kinetic energy), t is the time, p = pninnn 
is the mass density, p = 1.3 is the mean atomic mass 
(assuming cosmic abundances), tuh is the mass of the 
hydrogen atom, np is the hydrogen number density, Uq 
is the electron number density, u is the gas velocity, T is 
the temperature, and A(T) represents the radiative losses 
per unit emission measure. We used the ideal gas law, 
P = (7 — l)pe, where 7 = 5/3 is the adiabatic index. 

The calculations were performed using the flash 
code (jFrvxell et al.l l2QQQ[) . The hydrodynamic equa¬ 
tions were solved using the flash implementation 
of the Piecewise Parabolic Method (PPM) algorithm 
(|Colella fc WoodwardI Il984j) . For the present applica¬ 
tion, the code has been extended by additional computa¬ 
tional modules to handle the radiative losses A (through 
a table lookup/interpolation method) and to calcu¬ 
late the deviations from electron-proton temperature- 
equilibration and the deviations from equilibrium of ion¬ 
ization of the most abundant ions. The output of the 
latter modules was used in the synthesis of X-ray emis¬ 
sion (see Section E^l). 

The 3D simulations start once the majority of the en¬ 
ergy released in the explosion was kinetic (^ 26 hours 
after the explosion). Note that, at variance with our ID 
SN model, the 3D simulations do not include a heating 
term due to decays of radioactive isotopes synthesized in 
the SN explosion, such as ^^Co or ^"^Ti. On the other 
hand, in Appendix El we show that our simplification 
does not affect significantly the radial profiles of density 
and velocity of the ejecta at later times. This evidence 
supports our assumption to neglect the effect of radioac¬ 
tive heating during the interaction of the remnant with 
the H H region and with the dense equatorial ring. 

We followed the interaction of the blast wave and ejecta 
with the circumstellar nebula during the first 40 years of 
evolution (namely until th e presumed launch d ate of the 
Athena X-ray observatory: iNandr a et al.ll2Q13[) . The ini¬ 
tial remnant radius was ^ 20 AU (^ lQ~f. V>c). As sug¬ 
gested by previous studies (jOrlando et al.l[ 2012 f ). we as¬ 
sumed that the initial density structure of the ejecta was 
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Fig. 1.— A rendition in log scale of the circumstellar nebula 
around SN 1987A model initial conditions. The ring consisting of 
a uniform smooth component and high-density spherical clumps is 
shown in red; the HII region around the ring is marked by the blue 
clipped component. The white dot at the center of the coordinate 
system shows the position of the SN explosion. 


clumpy. To this end, after the ID radial density distribu¬ 
tion of ejecta was mapped into 3D, we modeled the ejecta 
clumps as per-cell random density pertur bations derived 
from a power-law probability distribution (|Orlando et al.l 
|2Q12[ ). The latter distribution, with index ^ = —1, is 
characterized by a parameter z^max representing the maxi¬ 
mum density perturbation allowed in the simulation. For 
the purposes of this paper, we assumed that the ejecta 
clumps have initial size 6 x 10^^ cm (0.4 AU, correspond¬ 
ing to 2% of the initial remnant radius) and a maximum 
density perturbation z/max = 5. The size and density 
contrast of the modeled ejecta clumps are in agreement 
with those suggested by spectropolarimetric studies of 
SNe (|Wang et al.ll2QQ3l 12004 [Hole et al']l2QlQ[ ). 

Initially the blast wave from the SN explosion 
propagates through the wind of the progenitor BSG. 
We adopted wind values disc ussed in the literature 
(|Morris fc Podsiadlowskil 120071 ): in particular, we as¬ 
sumed a spherically symmetric wind with a mass-loss 
rate of year“^ and wind velocity Uw = 

500 km s“^; the wind gas density is proportional to 
(where r is the radial distance from SN 1987A) and the 
termination shock of the wind is located approximately 
at Tw = 0.05 pc. Table [T] summarizes the parameters 
adopted. 

The circumstellar nebula probably originates from 
the interaction of a slow wind from an early RSG 
phase with a faste r wind from a subsequent BSG phase 
(jKwok et al.lll978l ). Thus the nebula was modeled as¬ 
suming that it is composed by a dense inhomogeneous 


equatorial ring surrounded by an ionized HIT regio n 
(|Chevalier fc DwarkadasI I1995 e ISugerman et al.l l2005f ). 
Figure [T] shows an example of initial configuration of the 
nebula. The ring (marked red in the figure) consists of 
a uniform smooth component and high-density spherical 
clumps mostly located in the inner portion of the ring. 
The smooth component has an elliptical cross section 
with the major axis w^g lying on the equatorial plane. 
The clumps mimic the protrusions emanating from the 
equatorial ring and probably formed by hydrodynamic 
instabilities caused by the interaction of a BSG and a 


RSG wind (|Sugerman et al.ir2QQ5f ). The clumps have all 
the same diameter Wc\ but their plasma density and ra¬ 
dial distance from SN 1987A are randomly distributed 
around the values < rid > and < Vd > respectively. The 
3D shape and geometry of the HII region were assumed 
to be analogous to those suggested from the analysis of 
HST data of the ring nebula around the B SG SBWl, a 
cand idate twin of the SN 1987A progenitor (|Smith et al.l 

i2fm . 

We explored the parameter space defined by: a) the 
plasma density nnii and inner edge thii of the HII re¬ 
gion; b) the density rirg, radial distance from SN 1987A 
Trg, radial extension rcrg, and height h^g of the uniform 
component of the ring; c) the average density < rid >, 
average radial distance from SN 1987A < rd >, diameter 
rCci, and number Nd of high-density spherical clumps of 
the ring. In order to limit as much as possible our explo¬ 
ration (that is very computer demanding in the case of a 
3D model), we adopted as fiducial values those reported 
in the literature and found by comparing the results of a 
ID hydrodynamic m odel of SN 1987A with observations 
(|Dewev et al.l l2Ql^ . Then we explored the parameter 
space around these values, adopting an iterative process 
of trial and error to converge on model parameters that 
reproduce the X-ray lightcurve and spectra of SN 1987A. 
Table [T] summarizes the ranges of values explored and 
the parameters of the model best-reproducing the data. 

We traced the evolution of the different plasma com¬ 
ponents (namely the ejecta, the HII region, and the ring 
material) by considering passive tracers associated with 
each of them. Each material is initialized with Ci = 1, 
while (7i = 0 elsewhere, where the index “i” refers to the 
ejecta (ej), the HII region (HII), and the ring material 
(rg). During the remnant evolution, the different mate¬ 
rials mix together, leading to regions with 0 < Ci < 1. 
At any time t the density of a specific material in a fluid 
cell is given by pi = pCi. 

We assumed the SN explosion at the origin of the 3D 
Cartesian coordinate system (xo,^0:^o) = (0^0? 0)- The 
system is oriented in such a way that the dense equa¬ 
torial ring lies on the (x, y) plane. The computational 
domain extends 1 pc in the x, and ^ directions. A ma¬ 
jor challenge in modeling the explosion and subsequent 
evolution of SN 1987A was the very small scale of the 
system in the immediate aftermath of the SN explosion 
(the initial remnant radius is ~ 10“^ pc) with the size of 
the rapidly expanding blast wave. To capture this range 
of scales, we exploited the adaptive mesh refinement ca¬ 
pabilities of the FLASH code, employing 18 nested levels 
of refinement, with resolution increasing twice at each 
refinement lev el. The refinement/derefinement criterion 
(iLohnerll 19871 ) implemented in flash follows the changes 
in mass density, temperature, tracer of ejecta, and tracer 
of the ring. 

The calculations were performed using also an auto¬ 
matic mesh derefinement scheme in the whole spatial 
domain that kept the computa tional cost approxirn ately 
constant as the blast expanded (jOrlando et al.ll2Q12[ ) : the 
maximum number of refinement levels used in the calcu¬ 
lation gradually decreased from 18 (initially) to 9 (at the 
final time) following the expansion of the blast and en¬ 
suring a number of grid zones per radius of the remnant 
>100 during the whole evolution. At the beginning (at 
the end) of the simulation, this grid configuration yielded 
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TABLE 1 

Adopted parameters of the CSM for the hydrodynamic models of the SN 1987A event. 


CSM component 

Parameters 

Units 

Range of values explored 

Best-fit values 

BSG wind: 

Mw 

{Mq year-1) 

10-7 

10-7 


Vw 

(km s-i) 

500 

500 


T w 

(pc) 

0.05 

0.05 

HII region: 

^HII 

(10^ cm“®) 

[0.8 - 3] 

0.9 


mil 

(pc) 

[0.08 - 0.2] 

0.08 

Equatorial ring: 

Urg 

(10^ cm-^) 

[1-2] 

1 


Trg 

(pc) 

0.18 

0.18 


Wrg 

(10^^ cm) 

[0.7 - 2] 

1.7 


hrg 

(IQl® cm) 

3.5 

3.5 

Clumps: 

< rici > 

(10"^ cm-^) 

[1-3] 

2.5 ±0.3 


< rci > 

(pc) 

[0.14 - 0.17] 

0.155 ±0.015 


'Wcl 

(10^® cm) 

[1-3] 

1.7 


Nci 


[40 - 70] 

50 


an effective resolution of 10“^ pc 5 x 10“^ pc) at 
the finest level, corresponding to ~ 100 zones per initial 
radius of the remnant (> 600 zones per final radius of 
the remnant). The effective mesh size varied from (10^)^ 
initially to 2048^ at the final time. Note that th e max - 
imum spatial resolution achieved by iPotter et al.l (|2014f ) 
in their 3D model was ^ Ax 10“^ pc (corresponding to an 
effective mesh size of 256^) during the whole evolution. 

The high-spatial resolution achieved in our simula¬ 
tions required about seven millions of CPU hours on 
the MareNostrum III cluster hosted at the Barcelona 
Supercomputing Center (Barcelona, Spain) and about 
four millions of CPU hours on the FERMI cluster hosted 
at CINECA (Bologna, Italy). Most of these resources 
were made available by a large computational program 
awarded in the framework of the PRACE^ (Partnership 
for Advanced Computing in Europe) initiative to enable 
high-impact scientific research with a pan-European su¬ 
percomputing infrastructure, the top level of the Euro¬ 
pean high performance computing systems. 

2.3. Synthesis of X-ray emission 

Erom the model results, we synthesized the X-ray emis¬ 
sion originating from the interaction of the blast wave 
with the surrounding nebula, followi ng an approach anal¬ 
ogou s to that of previous studies (jOrlando et al.l [20061 
l2QQ9f ). The results of numerical simulations are the evo¬ 
lution of temperature, density, and velocity of the plasma 
in the whole spatial domain. We rotated the system 
about the three axis to fit the inclinati on of the ring as 
found from the analysis of optical data (|Sugerman et al.l 
l2QQ5f) . namely U = 41°, iy = —8°, and iz = —9°. 

In the case of fast collisionless shocks, as those ob¬ 
served in SNRs, the synthesis of X-ray emission has to 
take into account the lack of temperature-equilibration 
between electrons and ions. In fact, the jump condi¬ 
tions at shock speeds exceeding 500 km s“^ drive ther¬ 
mal and dynamic changes of the plasma on very short 
times cales. Most of the sho ck energy is transferred to 
ions (|Ghavamian et aP 120071 ) , so that the ratio of the 
electron to ion temperature /3 in the post-shock region 
is generally less than I: the larger the shock velocity, 
the smaller /3. We computed this effect in our model 
by considering that the equilibration in the post-shock 


plasma is reached through Coulomb c ollisions which are 
expec ted to proceed relatively slowly (|Ghavamian et all 
leading to /3 < I. To take into account this effect, 
we added an additional passive tracer to the model equa¬ 
tions which stores the time tghj when the plasma in the 
j-th. domain cell was shocked. Then the electron temper¬ 
ature in each cell was calculated from the time tghj and 
from the ion temperature and plasma density derived by 
integrating Eqs. [H More specifically, we assumed that 
the electrons are heated almost istantaneously at the 
shock front up to kT ^0.3 keV (regar dless of the shock 
Mach number) by lower hybrid waves (|Ghavamian et ahl 
|2 oo 3), as suggested for shock velocities of the order of 
10^ km s“^ like those in our simulations. Then we cal¬ 
culated the evolution of ion and electron temperatures 
in each cell of the post-shock medium by including the 
effects of the Goulomb collisions. This evolution was cal¬ 
culated in the time Atj = t — tghj {t is the current time) 
since when the plasma in the cell was shocked. 

Another important effect to take into account in the 
synthesis of X-ray emission from fast shocks is the time 
lag of the plasma to change its ionization from a cool 
to a hot state. If the timescale of the temperature in¬ 
crease at the shock front is much shorter than the ion¬ 
ization and recombination timescales, the plasma ions 
can be at a lower ionization state than the equilibrium 
state corresponding t o the local electron temperature 
(jHamilton et al.l 1198311 . Such deviations from equilib¬ 
rium of ionization may have dramatic effects on the in¬ 
terpretation of observations. They are taken into ac¬ 
count in our mo del by following an app roach discussed 
in the literature (|Dwarkadas et alll2QIQ[ ) which is partic¬ 
ularly effective in the case of 3D simulations in order to 
have high efficiency together with a reasonable accuracy 
in the synthesis of X-ray emission. Erom the values of 
emission measure, electron temperature, and maximum 
ionization age in the j-th. domain cell, the corresponding 
X-ray spectrum is synthesized by using the NEI emission 
model VPSHOGK available in the XSPEG package along 
with the ”NEI vers ion 2.0” atomic data from ATOMDB 
(jSmith et al.lDOOU ). The emission measure in the j-th 
domain cell is emj = (where riej is the particle den¬ 
sity in the cell, is the cell volume); the electron tem¬ 
perature in the cell Tej is computed by taking into ac¬ 
count the deviations from temperature-equilibration as 
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TABLE 2 

Adopted parameters in the synthesis of X-ray emission 


Parameter 

Value^ 

Reference 

D 

51.4 kpc 

Panama f 19991 

Nu 

2.35 X 1021 cm-2 

Park et al. f20061 

He 

2.57 

Zhekov et al. f20091 

C 

0.03 

Zhekov et al. f20091 

N 

0.56 

Zhekov et al. f20091 

O 

0.081 

Zhekov et al. f20091 

Ne 

0.29 

Zhekov et al. f20091 

Mg 

0.28 

Zhekov et al. f20091 

Si 

0.33 

Zhekov et al. f20091 

s 

0.30 

Zhekov et al. f20091 

Ar 

0.537 

Zhekov et al. f20091 

Ca 

0.03 

Zhekov et al. f20091 

Fe 

0.19 

Zhekov et al. f20091 

Ni 

0.07 

Zhekov et al. f20091 


^ The abundances are rel ative to the solar photospheric values 
(| Anders &: GrevesselllQ^L 

described above; the maximum ionization age in the cell 
is Tj = n-ejAtj (where Atj is the time since when the 
plasma in the cell was shocked; see above). 

We assumed the source at a distance D = 51.4 kpc 
(|Panagialll999[ ). We adopted the metal abundances de¬ 
rived from the analysis of deep Chandra/LETG and 
HETG observations of SN 1987A (jZhekov et al.l l2QQ9[) 
and summarized in Table [2l The X-ray spectrum from 
each cell was filtered through the photoelectric absorp¬ 
tion by the ISM, assuming a col umn density A^h = 
2.35 X 10^^ cm“^ (jPark et al.ll2QQ6[ ). We integrated the 
absorbed X-ray spectra from the cells in the whole spa¬ 
tial domain. The spectra are then folded through the 
instrumental response of the X-ray instruments of inter¬ 
est, obtaining the relevant focal-plane spectra. In such 
a way, we derived X-ray spectra as they would be col¬ 
lected with XMM-Newton/EEIC and, in the near future, 
with Athena/WFl^ and X-ray images as they would be 
collected with Chandra/ACIS. The synthesized spectra 
and images were put in a format identical to that of real 
X-ray observations and analyzed with the standard data 
analysis system used for the specific instruments of in¬ 
terest. 

3. RESULTS 

3.1. Post-explosion evolution of the supernova 

The post-explosion evolution of the simulated 
SN I987A follows the t rend described in details by 
iPumo fc Zampieril (j2QIlf ). Here we summarize the main 
phases. The evolution is guided by the thermodynamics 
of the expanding ejecta and passes through three differ¬ 
ent phases. Initially the envelope is completely ionized 
and optically thick, and the emission is due to the re¬ 
lease of internal energy on a diffusion timescale (diffu¬ 
sive phase). Then the ejecta recombine and the emission 
is dominated by the sudden release of energy caused by 
the receding motion of the wavefront through the enve¬ 
lope (recombination phase). Finally, the envelope is re¬ 
combined and optically thin to optical photons, and the 
emission comes from the thermalization of the energy de¬ 
posited by gamma-ray photons (radioactive-decay phase 
or nebular phase). Usually the first two phases are glob- 


(a) 





Eig. 2.— Bolometric lightcurves (a) and evolution of the photo- 
spheric temperature (b) and velocity (c) derived from the models 
listed in Table [3] (solid c olored lines) comp ared to the corresponding 
quan t ities of SN 1987A (ICatchpole et aLlll987l . 119881 : IHamuv et ahl 

119881 : IPhillips et al.11198811 fempty squares and filled and empty 
diamonds). The photospheric temperature and velocity can be de¬ 
fined from the model during the so-called photospheric phase (see 
text), corresponding to the first ~ 100 days of evolution. Models 
SN-M17-E1.2-N8 and SN-M17-E1.2-N9 overlap each other. 

ally referred as photospheric phase, during which it is 
possible to compare the observed photospheric tempera¬ 
ture and velocity with the corresponding simulated quan¬ 
tities, contrarily to the nebular phase where the same no¬ 
tion of photosphere loses its meaning (|Pnmo fc Zampi^ 
l2QIll: l^rsten et al]l2QIl[) . We adopted the so-called in¬ 
flection time tinf (measured from the explosion epoch) as 
a measurement of the dur ation of the photospheric phase 
(|Pumo fc Zampieril |2Q 1 3[) : for SN I987A tinf ~ 100 days. 

From the models, we derived the main observables 
(namely the bolometric lightcurve and the evolution of 
the photospheric temperature and velocity) during the 
first 250 days after explosion. Then we compared the 
model results with observations of SN I987A by perform¬ 
ing a simultaneous fit of these observables, thus con¬ 
straining the ejected mass and the explosion energy of the 
SN. We found the comparison with the observations sat¬ 
isfactory for models with a total energy ranging between 
1.2 and 1.4 x 10^^ erg and an envelope mass between 15 
and 17 Mq (see Fig. [2]). Figure [3] shows the radial dis¬ 
tribution of mass density 26 hours after the breakout of 
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Fig. 3.— (Top) Radial density profiles of ejecta for the models 
reproducing the main observables of SN 1987A (see Table [3] and 
Fig. [2]) at t ~ 26 hours after the shock breakout. Different colors 
mark different models. (Bottom) Enlargment of the region marked 
with a box in the top panel, showing the details of the bulk of the 
ejecta. 


TABLE 3 

Parameters of the SN models reproducing the main 
OBSERVABLES OF SN 1987A DURING THE FIRST 250 DAYS FROM 
THE OUTBURST. 


Model 

■^env 

[Me] 

^SN 

[10^^ erg] 

a 

Rq 

[1012 cjjj] 

Mpfi 

[Me] 

SN-M17-E1.2-N8 

17 

1.2 

-8 

3 

0.07 

SN-M15-E1.2-N8 

15 

1.2 

-8 

3 

0.07 

SN-M17-E1.4-N8 

17 

1.4 

-8 

3 

0.07 

SN-M17-E1.2-N9 

17 

1.2 

-9 

3 

0.07 


the shock wave at the stellar surface for the models best 
reproducing the observations (see Table [3] for details). 
Note that runs SN-M17-E1.2-N8 and SN-M17-E1.2-N9 
differ only for the slope of the power-law describing the 
high-velocity shell of the SN (see Section [TT|) : the den¬ 
sity and kinetic energy in the outer envelope are lower 
for the steeper profile of density. The ejecta distribu¬ 
tion in this outer shell does not change significantly the 
observables of the SN. 

The difficulty of reproducing the observables at early 
time in Eig. [2] is mainly due to the initial conditions 
used in our simulations (see Section [2T|) . while some less 
prominent discrepancies at late times (80 — 110 days) 
might be explained by the absence of non-thermal ioniza- 
tion from gamma rays in our model (jPnmo fc Zampi^ 


HMH). 

The comparison of model results with observations en¬ 
abled us to constrain the bulk of envelope mass Menv and 
the enegy of the explosion E’sn- On the other hand, no 
firm conclusions on the distribution of energy and mass 
in the high-velocity shell of the SN were obtained from 
the analysis of SN observables. 

3.2. Remnant expansion through the eireumstellar 
nebula 

We followed the evolution from the SN phase to the 
SNR phase for 40 years, restricting our analysis to SN 
models reproducing the main observables of SN 1987A 
during the first 250 days of evolution (listed in Table [3]). 
Then, we explored the space of parameters describing the 
CSM (mainly the equatorial ring and the HII region) and 
compared the X-ray lightcurves and spectra synthesized 
from the models with those observed in SN 1987A (see 
Eig [5]). 

In all the cases investigated, the evolution follows the 
same general trend. Initially the blast wave from the 
SN propagates through the wind of the progenitor BSG. 
About 3 years after explosion, the ejecta reach the HII 
region (see on-line movie) and the transition from SN to 
SNR enters into its first phase of evolution (H Il-region- 
dominated phase). A forward and a reverse shocks are 
generated, the former propagating into the HII region 
and the latter driven back into the ejecta (see Eig. 11. 
The X-ray emission begins to increase rapidly with time 
and is dominated by emission from shocked plasma of 
the HII region and by a smaller contribution of shocked 
ejecta in the outer SN envelope. This is evident from 
the lightcurves (Eig. [5j), the spectra (Eig. [6]), and the 
emission morphology (Eig. [8] and on-line movie). The 
emitting plasma is largely out of equilibrium of ioniza¬ 
tion and its emission measure distribution as a function 
of electron temperature, /cTe, and ionization time scale, 
r, is peaked at kT^ « 2 keV and r ^ 3 x 10^^ s cm“^ 
with a sharp distribution that can be approximated with 
an isothermal plasma component (see upper panels of 
Eig. [1. This is in excellent agreement with the best-fit 
parameters derived from the spectral fitting with isother¬ 
mal components of X-ray spectra of SN 1987A at this 
epoch (see Appendix HI. This phase lasts until the blast 
wave hits the equatorial ring (year 15). 

In the H Il-region-dominated phase, we found that 
the sudden increase observed in the soft X-ray band 
([0.5, 2] keV) is best-fitted by models in which the outer 
layers of ejecta have a post-explosion radial profile of den¬ 
sity approximated by a power law with index a = —S (see 
Eig. [3)). On the contrary, we found that models with in¬ 
dex a = —9 systematically underestimate the soft X-ray 
flux during the first 7 years even by changing the values of 
nnii and rnii within the range of values compatible with 
observations. This is evident in Eig. [9] where we report 
the result for the model with a = —9 best approximating 
the observations: the discrepancy between model results 
and observations is evident for t < 7 yr even though the 
same model fits well the bolometric lightcurve of the SN 
during the first 250 days of evolution (see Eig. [a and the 
X-ray lightcurves for t > 7 years (Eig. [9]). Table [4] reports 
the parameters of this model to be compared with those 
of the model with a = —8 shown in Tabled! The reason 
for this discrepancy is the content of energy and mass of 
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Log n [ cm'^ ] 

2.0 2.5 3.0 3.5 4.0 ^ 


I. 

s? 


0.25 pc 
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t = 40 yr 

[0.5 - 2] keV 

2/2027 


% 


1 arcsec 



Fig. 4. — Interaction of the blast wave with the nebula. (Top) Three-dimensional volume rendering of particle density of the shocked 
plasma at the labeled times. (Middle) Corresponding synthetic maps of X-ray emission in the [0.5, 2] keV band integrated along the line-of- 
sight. Each image has been normalized to i ts maximum for vis ibility and convolved with a Gaussian of size 0.15 arcsec to approximate the 
spatial resolution of Chandra observations (|Helder et al.l[2M^ . (Bottom) Maps of X-ray emission of SN 1987A collected with Chandra at 
the labeled times, and normalized to their maximum for visibility (see Appendix I bI) . The overplotted ellipsoids represent the projection of 
circles lying in the equatorial plane of SN 1987A and fitting the position of the maximum X-ray emission in each observation. The dashed 
lines show an uncertainty of 10%. (See on-line movie). 


the outer shell of ejecta which is the first material hit¬ 
ting the nebula and, thus, determining the early X-ray 
emission from the ejecta-nebula interaction. Models with 
a = —9 have less mass and energy in the outer shell than 
models with a = —S and, as a result, the corresponding 
soft X-ray lightcurve rises more slowly and fails in fitting 
the observations. 

Note that the value of a have little (if any) effect in 
determining the main observables of the SN (e.g. the 
bolometric lightcurve) during the first 250 days of evolu¬ 
tion. SN models with the same envelope mass and total 
energy but with different values of a lead to very simi¬ 
lar results (see Fig. [2]). This is due to the fact that the 
observables of the SN depend on the bulk of ejecta. Our 
findings, therefore, shows that the X-ray emission orig¬ 
inating from the SNR in this phase can constrain the 
structure of outermost ejecta better than the emission 


from the SN. Since the density profile of ejecta is ex¬ 
pected to depend on the structure of the progenitor star 
and on the shock acceleration of the gas during the explo¬ 
sion, studying the early interaction of the ejecta with the 
nebula may provide important clues to the latest stage 
of stellar evolution and may be used as a probe of the 
mechanisms involved in the SN engine. 

During the H Il-region-dominated phase, the model en¬ 
abled us to constrain also the parameters characterizing 
the Hll-region (namely its density and inner radius); the 
best-fit parameters are listed in Table [H 

Around year 13 the blast wave hits the first dense 
clump of the equatorial ring and around year 15 it reaches 
the inner edge of the ring (see Fig. 0] and on-line movie). 
The transition from SN to SNR enters into the second 
phase of evolution (ring-dominated phase) which lasts 
until year 32, when the forward shock propagates be- 
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(a) Bolometric 





Years since SN 1987A 

Fig. 5.— Observed and modeled lightcurves. (a) Bolo¬ 
metric lightcurve of our favored model (cyan line) compared 

to the lightcnr ve _ of SN 19 8 7A (filled and em pty diamonds; 

ICatchpole et al ] I1987I . I1988I : IHamnv et al.l Il988lb (b) X-ray 
lightcurve in the [0.5,2] keV band synthesized from the favored 
model (black line) compared to the lightcurve of SN 19 87A ob¬ 
served with Rosat (grey diamonds; IHaberl e tjjjj2 006lb A SCA 
(brown; see Appendix IbI) . Chandra (magenta: IHgldgr_eL^ I^Qjjj) 
and XMM-Newton (cyan; IHaberl et al.ll200a : IMaggi et al.l 120121 : 
see Appendix IbI) . Green, blue and red lines mark the contribution 
to emission from the shocked ejecta, the shocked plasma from the 
HII region, and the shocked plasma from the ring, respectively, (c) 
Same as in Fig. [5h but for the lightcurve in the [3,10] keV band. 

yond the majority of the dense ring material. This phase 
is characterized by the interaction of the forward shock 
with the dense clumps of the ring. Each shocked clump 
evolves toward a core-plume structure with a crescent¬ 
like shape characterized by Kelvin-Helmholtz instabili¬ 
ties developing in the downstream region (see Fig.flOland 
on-line movie). As the shock travels through the clump, 
Rayleigh-Taylor instabilities develop on the upstream 
side of the clump, leading to its progressive fragmen¬ 
tation. A complex pattern of filaments and fragments 
forms in the interclump region, with densities varying 
between ~ 10^ cm“^ in the smooth component of the 
ring and ^ 10^ cm“^ in proximity of the clumps. Note 
that the appropriate description of the shock-clump in¬ 
teraction required the high-spatial resolution adopted in 
our simulations. The modelled features cannot be repro- 
duced at lower resolutions (e.g. iPotter et al.ll2Ql3 ). 


TABLE 4 

Parameters of the CSM for the hydrodynamic model with 

a = —Q BEST APPROXIMATING THE X-RAY LIGHTCURVES 


CSM component 

Parameters 

Units 

Best-fit values 

BSG wind: 

Mw 

(Mq year“i) 

10-'^ 


Vw 

(km s“^) 

500 


T "w 

(pc) 

0.05 

HII region: 

^HII 

(10^ cm“^) 

20 


rmi 

(pc) 

O.I 

Equatorial ring: 

Urg 

(10^ cm-^) 

I 


rrg 

(pc) 

0.18 


Wrg 

cm) 

1.7 


hrg 

cm) 

3.5 

Clumps: 

< rici > 

(ToWnU^) 

1.3 ±0.3 


< rd > 

(pc) 

0.17 ±0.015 


Wcl 

(IQl® cm) 

1.7 


Yi 


40 


In this phase, the contribution from the shocked ring 
dominates the X-ray emission (see Figs. [5]1H] and on-line 
movie). The shocked cores of the clumps lead predomi¬ 
nantly to soft X-rays and determine the further steepen¬ 
ing of the soft X-ray lightcurve (Fig. [Sb)- This emitting 
plasma component is essentially in collisional ionization 
equilibrium and its emission measure peaks to electron 
temperature kT^ « 0.5 keV and ionization time scale 
r 5 X 10^^ s cm“^ (see middle panels of Fig. 0. We can 
roughly identify this material with the plasma compo¬ 
nent with r > 10^^ s cm“^ derived from the spectral fit¬ 
ting of X-ray spectra of SN 1987A col lected with cmrent 
X-ray observatories for t > 15 yr (e.g. iHelder et al.ll2Q13l 
see also Appendix [B]) . The smooth component of the ring 
and the fragments of the shocked clumps stripped by hy¬ 
drodynamic instabilities dominante the emission in the 
hard band. This plasma component causes the broaden¬ 
ing of the emission measure distribution around the peak 
due to shocked clumps and includes plasma with kT^ up 
to ^ 2 keV and r down to ^ 10^^ s cm“^ (see middle 
panels of Fig. ED. Also a significant contribution to the 
emission with kT^ > 1 keV and r < 10^^ s cm“^ comes 
from the shocked ejecta. Given the complexity of the 
emission measure distribution in this phase, it is not ob¬ 
vious to attribute a physical meaning to the isothermal 
components with r < 10^^ s cm“^ derived from the spec¬ 
tral fitting of the X-ray spectra of SN 1987A and largely 
used in the literature (see Appendix [Bj. Indeed the phys¬ 
ical origin of the observed X-ray spectra is unveiled by 
our hydrodynamic model, as shown in Fig. [71 From the 
model best reproducing the X-ray lightcurves and spec¬ 
tra, we were able to constrain the parameters character¬ 
izing the equatorial ring (see Table [1] for details). 

During the first two phases of evolution, the rem¬ 
nant morphology in the X-ray band synthesized from the 
model appears very similar to that observed. In particu¬ 
lar, in the ring-dominated phase, the morphology is char¬ 
acterized by bright knots originating from the shocked 
clumps and resambles that of SN 1987A (see Figs. jH [51 
and on-line movie). In Fig. |4]we compare the synthetic 
and the observed morphology at two different epochs 
(years 14 and 26). The ellipsoids overplotted in the figure 
represent the projection of circles lying in the equatorial 
plane of SN 1987A and fitting the position of the maxi- 
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(c) 
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Fig. 6. — Synthetic and observed X-ray spectra of SN 1987A. (a) 
XMM-Newt on/EPlC-pn spectra at t = 14 yr. The true spectrum is 
marked in black (see Appendix I bI): the synthetic spectrum from the 
whole shocked plasma is marked in magenta; the contributions to 
emission from the different shocked plasma components are marked 
in green (ejecta), red (ring), and blue (HII region), (b) As in 
Fig. [6^, for t = 26 yr. (c) As in Fig. [GJi, for t = 40 yr and for 
spectra as they would be collected with Athena/WFl. 


mum X-ray emission in the observations. Although the 
extension of the X-ray source synthesized from the model 
seems to be smaller than the observed one (suggesting a 
modelled blast wave slightly slower than observed), the 
synthetic maps fit those observed within an uncertainty 
of 10% (dashed lines). In particular, at year 14, the ob¬ 
servations show a bright knot at north-west that may 
indicate that the observed blast wave is at a distance 


larger than that in our model. On the other hand, the 
knot is also well beyond the ellipsoid fitting the position 
of the forward shock in the equtorial plane, suggesting 
that, probably, the knot is the result of the interaction of 
the blast wave with some inhomogeneity at some height 
above the equatorial plane. This feature could be repro¬ 
duced in our model considering, for instance, an over- 
dense clump located well above the equatorial plane. 

The remnant enters into the third phase (ejecta- 
dominated phase) around year 32, when the contribution 
of shocked ejecta to the soft X-ray emission becomes the 
dominant component (see Figs. [MHl). The reverse shock 
travels through innermost ejecta with higher densities. 
Now the emitting plasma is characterized by a broad 
emission measure distribution that peaks at kT^ ^ 1 keV 
and T ^ 5x 10^^ s cm“^ although it is also characterized 
by few spikes with r > 10 ^^ s cm“^ due to the interaction 
of high-density clumps of ejecta with the ring (see lower 
panels of Fig. Eland right panel in Fig. [ 8 ]). The remnant 
morphology shows a revival of very bright knots but now 
due to shocked clumps of ejecta rather than clumps of the 
ring (see Fig.|4]and on-line movie). Our simulations show 
that, in this phase, the evolution of the X-ray emission 
depends again on the density distribution of the outer¬ 
most ejecta as in the H Il-region-dominated phase. Mod¬ 
els with the same explosion energy and envelope mass but 
with a different slope of the density profile of ejecta lead 
to significantly different X-ray lightcurves after year 32 
(compare Fig. [5] and Fig. [9]). In particular we found that 
the steeper is the density profile, the higher is the con¬ 
tribution of shocked ejecta to X-ray emission after year 
32. In the next future, therefore, SN 1987A will offer the 
possibility to study directly the structure and chemical 
composition of innermost ejecta and the imprint of the 
metal-rich layers inside the progenitor star. 


3.3. The central compact remnant 

Pioneering studies by iFransson fc Chevalierl (jl987[ ) 
suggested that the dense shells of expanding ejecta may 
obscure the X-ray emission of the compact object (if any) 
of SN 1987A for a few decades. Our model describes 
the distribution of the ejecta (shocked and unshocked) as 
well as the complex circumstellar environment. Thus the 
model allows us to get a thorough description of the local 
absorption in the remnant and to derive constraints on 
the X-ray emission of its yet undetected central source. 

The effective surface temperature of an isolated neu¬ 
tron star with an age of 10 — 30 y r is expected to be in 
the r ange Teff ^3 — 6 MK (e. g., IShternin fc YakovlevI 
l2008l ). This temperature should stay constant until the 
age of the star is lower than its thermal relaxation time, 
trei, which is in the range 30 — 300 yr, depending on 
the properties of the stell ar crust and of its core (e.g. 
IShternin fc YakovlevI 12008 ). From the analysis of Chan¬ 
dra observations. iNg et al.l (|2009t ) derived an upper limit 
of 0.010 counts per second for the central source (in 
the [0.08,10] keV energy band) on the basis of its non¬ 
detection in the Chandra-BRC data. Without account¬ 
ing for the local absorption, this upper limit corresponds 
to a surface temperature lower than ^2.5 MK, which is 
at odds with predictions from standard theories and may 
suggest extremely short values of trei due, for example, 
to superffuidity or quark matter. 

Our model shows that, at the present epoch, the col- 
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Fig. 7. — Left: distributions of emission measure vs. electron temperature kT^ and ionization time scale r at the labeled times, corre¬ 
sponding to the X-ray spectra shown in Fig. [6] Right: corresponding three-color composite images of the emission measure distributions. 
The colors show the contribution to emission measure from the different shocked plasma components, namely the ejecta (green), the ring 
(red), and the HII region (blue). 


umn density toward the center of the remnant^ is A^h ~ 
4.5—5 X10^^ cm“^, mainly associated with the unshocked 
ejecta. This local absorbing column density is almost a 

^ We inspected a region 3 x 10^^ cm wide (corresponding to 
0.2 arcsec considering the distance of 51.4 pc) to account for a 
possible proper motion of the central object around the geometrical 
center of the remnant, finding limited variations of the absorbing 
column. 


factor of 20 higher than the interstellar absorption to¬ 
ward the Large Magellanic Cloud and allows us to revise 
the constraints on t he therma l emiss ion from the central 
source suggested by iNg et all (|2QQ9f ) . By accounting for 
the local absorption, we then obtained an upper limit of 
3.9 MK for the surface temperature (i.e. a bolometric 
luminosity Lboi = 1.6 x 10^^ erg s“^), by considering a 
neutron star with radius Rns = 10 km at a distance of 
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Fig. 8.— Three-color composite images of the X-ray emission in the [0.5,2] keV band integrated along the line-of-sight at the labeled 
times. Each image has been normalized to its maximum for visibility and smoothed with a Gaussian of size 0.025 arcsec. The colors in the 
composite show the contribution to emission from the different shocked plasma components, namely the ejecta (green), the ring (red), and 
the HII region (blue). (See the on-line movie). 


(a) [0.5, 2.0] keV 




Fig. 9.— As in Fig. [Sj) and[5j: for a model with envelope mass 
Afenv = 17 Mq, ejecta energy Eq-^ = 1.2 x 10^^ erg, and with the 
density profile of ejecta in the high-velocity shell approximated by 
a power-law with index a = —9 (run SN-M17-E1.2-N9 in Table 
The parameters of the CSM of this model are reported in Table0] 

51.4 kpc. This revised value is in good agreement with 
standard theories, which sug gest bolometric luminositie s 
Tboi > 5 X 10^"^ erg s“^ (e.g. Yakovlev fc Pethi'c3l2QQ4l ). 
without invoking the presence of very short values of trei- 

The X-ray emission from young pulsars is typically 
dominated by a nonthermal component. Therefore, we 
also considered the case of nonthermal radiation, by as¬ 
suming a characteristic power-law index T = 1.5. Tak¬ 
ing into account the high local absorption from the 
unshocked ejecta, we found that the count-rate upper 
limit obtained with Chandra converts to a luminosity 
Lnt ^ 6 X 10^^ erg s“^ in the [2,10] keV band, cor¬ 
responding to a flux of ^ 0.1 mCrab at 51.4 kpc (to 
be compyed w ith Lnt ^ 7 x 10^^ erg s“^ obtained by 
iNg et, a,ijl2nn^ . 

Finally we studied the variation of the local absorp¬ 
tion with time, finding a relatively fast decrease. In 


Log n [ cm * 

2.0 2.5 3.0 3.5 

t = 26 yr 

1 


• 


1 


1 

0.05 pc 


Fig. 10. — Three-dimensional volume rendering of particle den¬ 
sity of the shocked plasma at t = 26 yr. The line-of-sight is aligned 
with the 2 ; axis and the supernova is located to the left of this plot. 
The figure shows a close-up view of the interaction of the ejecta 
with the ring in the equatorial plane. The bright knots are shocked 
clumps of the ring. (See the on-line movie). 


particular, we verified that the expansion of the ejecta 
will reduce the column density by a factor of ^ 2 
(Yh ~ 2 X 10^^ cm“^) at the end of our simulation (cor¬ 
responding to 40 yr after the explosion), namely at the 
presumed launch da te of the Athena X-ray observatory 
(|Nandra et al.|[^Q13[) . 

It is worth noting that, according to our model, the 
bulk of the absorption originates in the ejecta. The high 
metallicity of the expanding ejecta strongly enhances 
their optical depth with respect to a medium with a so¬ 
lar chemical composition, the O rich ejecta heavily ab¬ 
sorbing the X-ray radiation below 1 keV, and the Si, S, 
and Fe-rich ejecta contributing at higher energies (see 
IWilms et al.ll2QQQl for details). We can then consider our 
estimates (obtained by assuming a solar composition) as 
conservative. A detailed description of the local absorp- 
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tion, however, would requires a detailed knowledge of the 
distribution of the chemical abundances in the ejecta, 
which is beyond the scope of this paper. 

4. SUMMARY AND CONCLUDING REMARKS 

We investigated how the morphology and the emission 
properties of the remnant of SN 1987A reflect: 1) the 
physical characteristics of the progenitor SN and 2) the 
early interaction of the SN blast wave with the surround¬ 
ing inhomogeneous nebula. To this end, we have devel¬ 
oped a model describing SN 1987A from the breakout of 
the shock wave at the stellar surface up to the 3D expan¬ 
sion of its remnant. A major challange was capturing the 
enormous range in spatial scales (spanning six orders of 
magnitude) that required a very high spatial resolution 
(down to ^ 0.2 AU). We performed an exploration of 
the parameters space, searching for a model reproduc¬ 
ing at the same time both the observables of the SN 
(i.e. bolometric lightcurve, evolution of line velocities, 
and continuum temperature at the photosphere) and the 
X-ray emission of the remnant (lightcurves, spectra, and 
morphology). Our findings lead to several conclusions: 

1. We identified three phases in the evolution (see FigH] 
and on-line movie). During the first phase (Hll-region- 
dominated phase; from year 3 to 15) the fastest ejecta 
interact with the HII region and the X-ray emission is 
dominated by shocked plasma from this region and by 
a smaller contribution from the outermost ejecta (see 
Figs.[5l[6l and on-line movie). In the second phase (ring- 
dominated phase; from year 15 to 32) the ejecta interacts 
with the dense equatorial ring. The emission in the soft 
X-ray band is largely dominated by the shocked clumps 
(see Figs. [5l [6l and on-line movie). The emission in the 
hard X-ray band is mainly due to shocked plasma from 
the smooth component of the ring and to fragments of 
the shocked clumps stripped by hydrodynamic instabil¬ 
ities developing at the cloud boundaries. In the third 
phase (ejecta-dominated phase; from year 32), the for¬ 
ward shock propagates beyond the majority of the dense 
ring material and the reverse shock travels through the 
inner envelope of the SN. The X-ray emission is domi¬ 
nated by shocked ejecta (see Figs. [5] and [6]) and the rem¬ 
nant morphology is characterized by very bright knots 
due to shocked clumps of ejecta (see on-line movie). 

2. Our favored model reproduces altogether the main 
observables of SN 1987A (bolometric lightcurve and evo¬ 
lution of photospheric temperature and velocity) during 
the first 250 days after explosion and the observed X- 
ray lightcurves and spectra of its remnant (see Figs. [5] 
and © in the following 30 years. Therefore, we have 
demonstrated that the physical model reproducing the 
observables of the SN is able to reproduce also the observ¬ 
ables of the subsequent expanding remnant, providing a 
firm link between two research fields (SN explosions and 
SNR evolution) which, traditionally, are based on models 
that are independent from each other. In other words, 
in the case of SN I987A, we have demonstrated the con¬ 
sistency between the cause (the SN explosion) and the 
effect (the interaction of the remnant with the surround¬ 
ing medium). This is a great advance with respect to 
a parameterised explosion model to bootstrap the SNR 
(as those commonly used in the literature) whose pa¬ 
rameters are chosen to fit the observables of the remnant 
but that does not ensure to fit also the observables of the 


progenitor SN. 

3. From our favored model, we identified the im¬ 
print of SN I987A on the remnant emission. In par¬ 
ticular, we constrained the SN explosion energy in the 
range 1.2 — 1.4 x 10^^ erg and the envelope mass in the 
range 15 — I7M0. The model also constrained the phys¬ 
ical properties of post-explosion ejecta. During the HII- 
region-dominated phase, the sudden increase of X-ray 
flux observed around year 4 (Fig. [Sb) is reproduced if 
the outermost ejecta have a post-explosion radial pro¬ 
file of density approximated by a power law with index 
a = —S. On the contrary, mod els with index a = — 9 (as 
early suggested for SN I987A; IWooslevlll988[ ) systemat¬ 
ically underestimate the soft X-ray flux during the first 
7 years since the explosion (see Fig. E), independently of 
the density structure of the nebula within the range of 
values compatible with observations. Indeed the shape 
of the lightcurves in this phase reflects the structure of 
outer ejecta and reveals the imprint of the SN on the 
remnant emission. 

4. Our favored model allowed us to constrain the 
structure of the pre-SN nebula. In the ring-dominated 
phase, the shape of lightcurves and spectra reflect the 
density structure of the nebula, allowing to disentan¬ 
gle the effects of the SN event (identified in the previ¬ 
ous phase) from those of remnant interaction with the 
environment. This enabled us to ascertain the origin 
of the multi-thermal X-ray emission and to constrain 
the nebula structure (see Table [T] for all the details). 
From this model we estimated that the total mass of 
the ring is Mrg = 0.062 Mq of which ^ 64% is plasma 
with density n ~ 1000 cm“^ and ^ 36% is plasma with 
n ~ 2.5 X 10^ cm“^. These values are in excellent agree¬ 
ment with those derived from optical spectroscopic data 
for the density stru cture of the ionized gas of the ring 
(|Mattila et al.l l l2ninD . 

Our model enabled us to make predictions about the 
ejecta evolution and the future changes in the remnant 
morphology in X-rays. In the next few years, the rem¬ 
nant will enter in the ejecta-dominated phase. The X- 
ray flux will reflect the radial profile of density in outer 
ejecta: the steeper the slope of this profile, the higher 
the emission from shocked ejecta. The emission will en¬ 
able us to study in more details the ejecta asymmetries 
and the distribution of metal-rich layers. This will pro¬ 
vide important clues on the dynamics of the explosion 
and might even improve our knowledge about the nu¬ 
cleosynthesis processes occurring in the latest stage of 
stellar evolution and during the SN explosion, making 
the remnant of SN I987A a unique probe of CC-SNe. 

Finally, we investigated how our model relate to the ex¬ 
istence of the yet undetected central compact remnant. 
The complete picture of the line of sight column towards 
the center of SN I987A provied by our model has shown 
that the emission of the compact remnant cannot be re¬ 
vealed yet due to a local absorbing column density which 
is a factor of 20 higher than the interstellar absorption to¬ 
ward the Large Magellanic Cloud. The constraint on the 
thermal emission from the central source inferred from 
our model is in good agreement with standard theories 
of neutron star cooling. 
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Fig. A1. — Radial profiles of mass density (upper panel) and velocity (lower panel) of ejecta at t = 1 yr since the SN event. The solid 
curves show the angle-averaged profiles derived with the 3D hydrodynamic model describing the SNR evolution; the dashed curves show 
the analogous profiles derived with the ID relativistic radiation hydrodynamic code describing the post-explosion evolution of the SN. 

APPENDIX 

EFFECT OF RADIOACTIVE HEATING DURING THE REMNANT EXPANSION 

The 3D hydrodynamic simulations describing the evolution of the SNR and its interaction with the inhomogeneous 
CSM do not include a heating term due to decays of radioactive isotopes synthesized in the SN explosion (e.g. ^^Co 
or ^"^Ti). To check the validity of our assumption, we used our ID relativistic radiation hydrodynamic code (which 
includes the radioactive heating term; see Sect. 1^ . to extend run SN-M17-E1.2-N8, simulating the post-explosion 
evolution of the SN, and covering the first year of evolution. Then we compared the radial profiles of mass density 
and velocity of ejecta derived in this case at t = 1 yr with the angle-aver aged radial profiles of density an d ve locity 
derived at the same epoch with our 3D hydrodynamic model (not including the radioactive heating). Fig. I All shows 
some differences in the density profiles that are mainly due to the effect of ejecta clumping considered only in the 
3D simulation. On the other hand, the two velocity profiles are very similar suggesting that the effect of radioactive 
heating does not affect too much the dynamics of the ejecta. This result suggests that we can safely neglect the effect 
of radioactive heating during the interaction of the remnant with the surrounding nebula. 

THE DATA ANALYSIS 

We present new results of the analysis of different observations of the remnant of SN 1987A performed with the 
ASCA^ XMM-Newton^ and Chandra X-ray telescopes. We adopted Xselect V2.4, CIAO V4.6.1, and SAS V13 for the 
reduction of ASCA/SIS, XMM-Newton/EFIC, and Chandra/ACIS data, respectively. All the spectra were analyzed 
with XSPEC V12.8.2. 


ASCA data 

We analyzed the ASCA observations ID 55039000 (performed on November 1997), 56041000 (November 1998), and 
57034000 (November 1999), to derive the X-ray unabsorbed fluxes in the [0.5,2] keV and [3,10] keV bands, shown in 
Fig. m We screened the data according to the standard screening criteria and we added the screened SISO and SISl 
spectra by using the ADDASCASPEC task. We verified that the spectra extracted from observations 55039000 and 
56041000 were consistent with each other, considering the relatively low signal-to-noise ratio of the early observations 
of SN 1987A. We then fitted the SIS spectra extracted from observations 55039000 and 56041000 simultaneously. For 
each spectrum, we subtracted a background spectrum extracted from a nearby region immediately outside of the 
supernova remnant (SNR) shell and we verified that the best-fit values do not depend significantly on the choice of 
the background region. Spectra were modeled by an absorbed (TBABS model in XSPEC) optically thin plasma in 
non-equilibrium of ionization (VNEI model). In the fittings, the plasma chemical abundances were fixed to the values 
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derived from the analysis of deep Chandra/ yETG and HET G observations of SN 1987A (|Zhekov et al.|[2QQ^ . while 
the interstellar column density was fixed to (jPark et alJl2QQ6[ ) A^h = 2.35 x 10^^ cm“^. Our best-fit models provide a 
good description of the spectra = 41.5 with 47 d.o.f. for the 1997 and 1998 spectra, and = 36.0 with 29 d.o.f. 
for the 1999 spectrum) and the best-fit temperatures and ionization timescales are kT = 1.5±0.3 keV, r = 3±2 x 10^^ 
s cm“^ for the 1997 and 1998 spectra and kT = 1.7 ± 0.4 keV, r = 2 ± 1 x 10^^ s cm“^ for the 1999 spectrum. 

XMM-Newton data 

As for the XMM-Newton data, we analyzed observation ID 0083250101 (performed on April 2001) and the previously 
unpublished observation ID 0690510101 (performed on December 2012). We focussed on the EPIC data and we 
selected events with PATTERN< 12 for the MOS cameras, PATTERN<4 for the pn camera, and ELAG=0 for both. 
We screened the original event files by using the sigma-clipping algorithm (ESPEILT tasks). The pn spectra extracted 
from the 2001/2012 observations are shown in the upper/middle panels of Eig.[6l To derive the [0.5, 2] keV and [3,10] 
keV fluxes of the 2012 observation, we fitted the pn and MOS spectra simultaneously (we selected a circular extraction 
region with radius r = 20") in the [0.3,9] keV energy band. We adopte d a model includin g three components of 
optically thin plasma in non-equilibrium of ionization (NEI), widely used (iMaggi et al.ll2012h t o describe the latest 
X-ray observations of SN 1987A. Again we fixed the interstellar column density (iPark et al.ll2006l) to /Vh = 2.35 x 10^^ 
cm“^. The chemical abundances were fixed to the values derived in previous studies (jzhekov et al.ll200^ for all the 
three components. The best-fit temperatures are kTi = 0.48lo!o2 ^^2 = 0.77 ± 0.03 keV, and kTs = 2.41 ± 0.05 

keV, while the corresponding ionization timescales are ri = 1.5 ± 0.1 x 10^^ s cm“^, r 2 = 2.4 ± 0.3 x 10^^ s cm“^, 
and rs = l.lltg;^^ X 10^^ s cm The presence of a component with a r > 10^^ s cm ^ (in our case, component 2) 
indicates that part of t he X-ray emitting p lasma has reached the collisional ionization equilibrium and is in agreement 
with previous findings (|Helder et al.ll201^ . We point out that these models only provide heuristic (phenomenological) 
descriptions of the spectra, whose physical origin can be unveiled only by accurate hydrodynamic modeling (see 

Section [321) • 


Chandra data 

Einally, we analyzed the Chandra/XCl^ observation 1967 (performed on December 2000) and the previously un¬ 
published observation 14697 (March 2013) to produce the X-ray images shown in Eig. |4] (lower panels). To study 
the morphology of SN 1987A from th e observations, we carefully followed the procedure described in the literature 
(jBurrows et al.lDOOOUPark et al.ll2002r). T he ACIS data were deconvolved with a maximum likelihood algorithm (with 
25 iterations: lRichardsonlll972l: lLucvlll974l ). using an on-axis point-spread function produced by the MARX simulation 
package. The high photon statistics allowed us to use 0.062" pixels for the deconvolution process. The deconvolved 
images were finally smoothed with a Gaussian with a = 0.1". 
















